#include "StdAfx.h"

#include "EphemDE.h"
#include "CustomException.h"


CBufferHolder::CBufferHolder(FILE * file, int nBufSize)
: m_File(file)
{
	if(m_File == NULL)
		throw CCustomException(_T("CBufferHolder::CBufferHolder: Empty File"));
	if(nBufSize == 0)
		throw CCustomException(_T("CBufferHolder::CBufferHolder: Zero size"));

	for(int i = 0; i < sizeof(m_buffer)/sizeof(m_buffer[0]); i++)
	{
		m_buffer[i].last = (i == 0)? true: false;
		m_buffer[i].Offset = -1;
		m_buffer[i].vec.resize(nBufSize);
	}
}

double * CBufferHolder::GetBufferAtOffset(int nOffset)
{
	if(nOffset < 0)
		throw CCustomException(_T("CBufferHolder::GetBufferAtOffset: Offset less then zero"));

	for(int i=0; i < sizeof(m_buffer)/sizeof(m_buffer[0]); i++)
	{
		if(	m_buffer[i].Offset == nOffset)
			return &(m_buffer[i].vec[0]);
	}

	Buffer * buf = (m_buffer[0].last)? (m_buffer+1): (m_buffer);
	fseek(m_File, nOffset, SEEK_SET);

	size_t i = fread(&(buf->vec[0]), sizeof(buf->vec[0]), buf->vec.size(), m_File);
	if (i < buf->vec.size() && feof(m_File)) 
	{
		throw CCustomException(_T("CBufferHolder::GetBufferAtOffset: error reading: read %d koef from %d%d"),
		i, buf->vec.size(), (feof(m_File))? " EOF": "");
	}
	else
	{
		for(int i=0; i < sizeof(m_buffer)/sizeof(m_buffer[0]); i++)
			m_buffer[i].last = false;
		buf->last = true;
		buf->Offset = nOffset;
		return &(buf->vec[0]);
	}
}

// ====================================== //

CEphemDE::CEphemDE(LPCTSTR stzEphemFileName)
: m_file(stzEphemFileName, _T("rb"))
, m_lastx(0.0)
, m_lastxv(0.0)
{
	ephcom_readbinary_header ();

	m_ReadBlockBuff = CBufferHolder::AP(new CBufferHolder(m_file, m_header.ncoeff));

	m_vecPos.resize(m_header.maxcheby);
	m_vecVel.resize(m_header.maxcheby);

	LPCSTR stzName [] = { 
		("GM1"), 
		("GM2"), 
		("GMB"),
		("GM4"), 
		("GM5"), 
		("GM6"), 
		("GM7"), 
		("GM8"), 
		("GM9"), 
		("EMRAT"), 
		("GMS") };

		if(sizeof(stzName)/sizeof(stzName[0]) != sizeof(m_dGM)/sizeof(m_dGM[0]))
			throw CCustomException(_T("wrong BODY_COUNT value, must be 12"));

		for(int i = 0; i<BODY_COUNT; i++)
		{
			m_dGM[i] = GetParam(stzName[i]);
		}

		m_dGM[eMOON] = m_dGM[eEARTH] / (1.0 + m_dGM[eMOON]);
		m_dGM[eEARTH] -= m_dGM[eMOON] ;
}

CEphemDE::~CEphemDE(void)
{
}

bool CEphemDE::GetCoord(double date, bool bBaryCenter, AllBodies & OutCoords, bool bGetVelocity)
{
	if(date <= m_header.ss[0])
		throw CCustomException(_T("date less then ephemerid range"));

	if(date > m_header.ss[1])
		throw CCustomException(_T("date great then ephemerid range"));

	memset( &OutCoords, 0, sizeof(AllBodies));

	struct ephcom_Coords coords;
	coords.bary = bBaryCenter? 1: 0;
	coords.center = EPHCOM_SUN;
	coords.km = 0;        /* Use AU, not kilometers */
	coords.seconds = 0;   /* Timescale is days, not seconds */
	coords.et2[0] = date;  /* Good enough precision for test dates */
	coords.et2[1] = 0.0;

	if ( ephcom_get_coords(&coords, bGetVelocity) ) 
		return false;

	for(int i = 0; i<BODY_COUNT; i++)
	{
		OutCoords.bodies[i].GM = m_dGM[i];
		memcpy(OutCoords.bodies[i].x, coords.pv[i]    , sizeof(OutCoords.bodies[i].x) );
		if(bGetVelocity)
			memcpy(OutCoords.bodies[i].v, coords.pv[i] + 3, sizeof(OutCoords.bodies[i].v) );
	}

	return true;
}

double CEphemDE::GetParam(LPCSTR szName)
{
	int len = strlen(szName);
	if(len > 6)
		throw CCustomException(_T("too long"));

	for(int pos = 0; pos < sizeof(m_header.cnam)/sizeof(m_header.cnam[0]); ++pos)
	{
		if(memcmp(m_header.cnam[pos], szName, len*sizeof(szName[0])) == 0)
			return m_header.cval[pos];
	}

	throw CCustomException(_T("can not found param %s"), szName);
}

// =================================================================== //

double ephcom_indouble(FILE *infp) 
{
	double retval;
	int wasread = fread(&retval, sizeof(double), 1, infp);
	return(retval);
}

int ephcom_inint(FILE *infp) 
{
	int retval;
	fread(&retval, sizeof(retval), 1, infp);
	return retval;
}

// =================================================================== //

/*
Read a JPL Ephemeris header in binary format.  Store values in
an ephcom_Header struct.
*/
int CEphemDE::ephcom_readbinary_header() 
{
	int i, j;
	FILE * infp = m_file;
	ephcom_Header *header = &m_header;

	if (infp != stdin) 
		rewind(infp);
	/*
	Read title lines.
	*/
	for (i=0; i<3; i++) 
	{
		for (j=0; j<84; j++) 
		{
			header->ttl[i][j] = fgetc(infp);
		}
		if (i == 0 && strncmp(header->ttl[0], "JPL ", 4) != 0) 
		{
			fprintf(stderr,"\nERROR: file is not a JPL ephemeris file.\n\n");
			if (strncmp(header->ttl[0], "KSIZE", 5) == 0)
				fprintf(stderr,"File is an ASCII JPL ephemeris header instead.\n\n");
			exit(1);
		}
		header->ttl[i][j] = '\0';
	}
	/*
	Read constant names.
	*/
	for (i=0; i<400; i++) 
	{
		for (j=0; j<6; j++) 
		{
			header->cnam[i][j] = fgetc(infp);
		}
		header->cnam[i][j] = '\0';
	}
	/*
	Read ephemeris start epoch, stop epoch, and step size (in Julian Days).
	*/
	for (i=0; i<3; i++) 
	{
		header->ss[i] = ephcom_indouble(infp);
	}
	/*
	Read NCON, AU, EMRAT.
	*/
	header->ncon  = ephcom_inint(infp);
	header->au    = ephcom_indouble(infp);
	header->emrat = ephcom_indouble(infp);
	header->nval  = header->ncon;
	/*
	Read indexes for coefficients in data block.  Written in transposed
	order (Fortran and C matrices are transposed).
	*/
	for (i=0; i<12; i++) 
	{
		for (j=0; j<3; j++) 
		{
			header->ipt[i][j] = ephcom_inint(infp);
		}
	}

	header->numde = ephcom_inint(infp);  /* Get ephemeris number */

	for (i=0; i<3; i++) 
		header->lpt[i] = ephcom_inint(infp);
	
	/*
	If there are no coefficients for an ipt[i][] object (i.e., ipt[i][1]==0),
	then ipt[i][0] should contain the value of the next available coefficient
	number rather than 0, as per communication of Myles Standish to Paul Hardy
	on preferred format of ephemeris headers.

	If there are no libration coefficients (i.e., lpt[1]==0), then lpt[0]
	should contain the value of the next available coefficient number rather
	than 0 as well, as per the same communication from Myles Standish.
	*/
	/* First set j to maximum index into ipt[] that has coefficients */
	j = 0;
	for (i=1; i<12; i++)
	{
		if (header->ipt[i][1] > 0 && header->ipt[i][0] > j)
			j = i;
	}

	/* Now set j to next available index count. */
	if (header->lpt[1] > 0 && header->lpt[0] > j)
		j = header->lpt[1] + header->lpt[1] * header->lpt[2] * 3;
	else
		j = header->ipt[j][0] +	header->ipt[j][1] * header->ipt[j][2] * (j==11 ? 2 : 3);

	for (i=1; i<12; i++)
	{
		if (header->ipt[i][0] == 0) 
			header->ipt[i][0] = j;
	}
	if (header->lpt[0] == 0)
		header->lpt[0] = j;

	/*
	Set the maximum number of Chebyshev coefficients possible for this file,
	to initialize position and velocity Chebyshev coefficient arrays during
	Chebyshev interpolation.
	*/
	header->maxcheby = 0;
	for (i=0; i<12; i++)
	{
		if (header->ipt[i][1] > header->maxcheby)
			header->maxcheby = header->ipt[i][1];
	}

	if (header->lpt[1] > header->maxcheby)
		header->maxcheby = header->lpt[1];

	// /*
	//    From JPL ephemeris number, set NCOEFF and calculate KSIZE = 2*NCOEFF.
	// */
	// switch (header->numde) {
	//    case 102:
	//       header->ncoeff = 773;
	//       break;
	//    case 200:
	//       header->ncoeff = 826;
	//       break;
	//    case 202:
	//       header->ncoeff = 826;
	//       break;
	//    case 403:
	//       header->ncoeff = 1018;
	//       break;
	//    case 405:
	//       header->ncoeff = 1018;
	//       break;
	//    case 406:
	//       header->ncoeff = 728;
	//       break;
	//    case 410:
	//       header->ncoeff = 1018;
	//       break;
	//    default:
	//       header->ncoeff = 1018;
	//       break;
	//    }
	/*
	Calculate number of coefficients, starting with
	highest index into a data block (stored in j).
	*/
	j = 0;
	for (i=1; i<12; i++)
	{
		if (header->ipt[i][1] > 0 && header->ipt[i][0] > header->ipt[j][0]) j = i;
	}
	/*
	Now see if the starting point we found is lower than where
	lpt[] starts.  If not, use lpt[] for largest value.
	*/
	if (header->lpt[1] > 0 && header->lpt[0] > header->ipt[j][0]) 
	{
		header->ncoeff = header->lpt[0] - 1 +     /* starting point */
			(header->lpt[1] *        /* coefficients per coordinate */
			header->lpt[2]) *       /* subblocks per block */
			3;                      /* coordinates */
	}
	else 
	{
		header->ncoeff = header->ipt[j][0] - 1 +  /* starting point */
			(header->ipt[j][1] *     /* coefficients per coordinate */
			header->ipt[j][2]) *    /* subblocks per block */
			(j == 11 ? 2 : 3);       /* coordinates */
	}
	// printf("Number of coefficients: %d\n", header->ncoeff);

	header->ksize = header->ncoeff + header->ncoeff; /* KSIZE = 2*NCOEFF always */
	/*
	Skip to second block in file.
	*/
	fseek(infp, header->ncoeff * 8, SEEK_SET);
	/*
	Read ephemeris constants.
	*/
	for (i=0; i<header->ncon; i++) {
		header->cval[i] = ephcom_indouble(infp);
		if (strncmp(header->cnam[i], "NCOEFF", 6) == 0) 
		{
			header->ncoeff = (int)header->cval[i];
			header->ksize  = (int) (2 * header->cval[i]);
			_ASSERT(header->ncoeff == header->cval[i]);
			_ASSERT(header->ksize  == (2 * header->cval[i]));
		}
		else if (strncmp(header->cnam[i], "LENUM ", 6) == 0)
			header->numle = (int) header->cval[i];
	}
	if (header->numle == 0) header->numle = header->numde;

	return(0);
}

/*
ephcom_get_coords() - Interpolate positions and velocities at given time.
*/
int CEphemDE::ephcom_get_coords(struct ephcom_Coords *coords,	bool bGetVelocity) 
{
	FILE * infp = m_file;
	ephcom_Header *header = &m_header;
	//double *datablock = &m_vecDataBlock[0];

	double et2[2];    /* Ephemeris time, as coarse (whole) and fine time  in JD */
	double totaltime; /* Sum of whole and fractional JD */
	double filetime;  /* JDs since start of ephemeris file */
	double blocktime; /* JDs since start of data block */
	double subtime;   /* JDs since start of subinterval in block */

	int i, j;
	int blocknum;
	int subinterval; /* Number of subinterval for this body */
	int dataoffset; /* Offset in datablock for current body and subinterval */
	double subspan; /* Span of one subinterval in days */
	double chebytime; /* Normalized Chebyshev time, in interval [-1,1]. */
	int ncoords; /* Number of coordinates for position and velocity */
	int retval; /* Return value */

	retval = 0; /* Assume normal return */
	/*
	Split time JD into whole JDs (et2[0]) and fractional JD (et2[1]).
	*/
	totaltime = coords->et2[0] + coords->et2[1];
	if (totaltime <= header->ss[0] || totaltime >= header->ss[1]) 
	{
		// fprintf(stderr,"Time is outside ephemeris range.\n");
		retval = -1;
	}
	else 
	{
		et2[0] = (int)totaltime;
		et2[1] = (coords->et2[0] - et2[0]) + coords->et2[1];
		filetime = totaltime - header->ss[0]; /* Days from start of file */
		blocknum = (int)(filetime / header->ss[2]); /* Data block in file, 0.. */
		/*
		Read the data block that contains coefficients for desired date
		*/
		long fileOffset = (blocknum + 2) * header->ncoeff * sizeof(double); /* 8 bytes per coefficient */
		double *datablock = m_ReadBlockBuff->GetBufferAtOffset(fileOffset);

		// printf("%12.2f : Block %d, %12.2f to %12.2f\n",
		//        totaltime, blocknum, datablock[0], datablock[1]);
		// ephcom_parse_block(stdout, header, datablock);
		/*
		Now step through the bodies and interpolate positions and velocities.
		*/
		blocktime = totaltime - datablock[0]; /* Days from block start */
		for (i=0; i<13; i++) 
		{
			//    printf("Calculating Body %d\n", i+1);
			subspan = header->ss[2] / header->ipt[i][2]; /* Days/subinterval */
			subinterval = (int)((totaltime - datablock[0]) / subspan);
			if (coords->seconds) subspan *= 86400.0; /* for km/sec, 86400 sec/day */
			//    printf("subspan = %g\n", subspan);
			ncoords = i == 11 ? 2 : 3; /* 2 coords for nutation, else 3 */
			//    printf("Body %d for JD %12.2f subinterval=%d [%d %d %d]\n",
			//           i, totaltime, subinterval,
			//              header->ipt[i][0], header->ipt[i][1], header->ipt[i][2]);
			dataoffset = header->ipt[i][0] - 1 +
				ncoords * header->ipt[i][1] * subinterval;
			//    printf("Body %d, Subinterval %d, Dataoffset = %d\n",
			//           i+1, subinterval, dataoffset);
			subtime = blocktime - subinterval * header->ss[2] / header->ipt[i][2];
			/*
			Divide days in this subblock by total days in subblock
			to get interval [0,1].  The right part of the expression
			will evaluate to a whole number: subinterval lengths are
			all integer multiples of days in a block (all powers of 2).
			*/
			chebytime = subtime / subspan;
			chebytime = (chebytime + chebytime) - 1.0;
			if (chebytime < -1.0 || chebytime > 1.0) {
				fprintf(stderr,"Chebyshev time is beyond [-1,1] interval.\n");
				fprintf(stderr,
					"filetime=%f, blocktime=%f, subtime=%f, chebytime=%f\n",
					filetime, blocktime, subtime, chebytime);
			}
			else {
				/*
				Everything is as expected.  Interpolate coefficients.
				*/
				ephcom_cheby(chebytime, subspan, &datablock[dataoffset], ncoords, m_header.ipt[i][1], coords->pv[i], bGetVelocity);
			}
		}
		// /*
		//    Set any user-defined coordinates to zero.
		// */
		// for (i = 16; i < EPHCOM_NUMOBJECTS; i++)
		//    coords->pv[i][0] = coords->pv[i][1] = coords->pv[i][1] =
		//       coords->pv[i][1] = coords->pv[i][1] = coords->pv[i][1] =  0.0;
		/*
		With interpolations complete, calculate Earth from EMBary and
		Sun from SSBary.  Preserve other coordinates.
		*/

		// 		for (j=0; j<6; j++) 
		// 			printf("%25.20lf\n", coords->pv[1][j]/header->au);
		// 		printf("\n");
		// 		for (j=0; j<6; j++) 
		// 			printf("%25.20lf\n", coords->pv[2][j]/header->au);

		if (!coords->km)  /* Calculate AU, not kilometers */
		{
			double au1 = 1.0/header->au;
			for (i=0; i<15; i++) 
			{
				double *PVal = &(coords->pv[i][0]);
				if (i == 13) 
					i = 15; /* Skip over nutations and librations */
				for (j=0; j<6; j++, PVal++ )
				{
					*(PVal) *= au1;
					//*(PVal) /= header->au;
				}
			}
		}

		for (j=0; j<6; j++) 
		{
			coords->pv[15][j] = coords->pv[ 9][j]; /* Save original lunar coords */
			coords->pv[14][j] = coords->pv[12][j]; /* Librations if on file */
			coords->pv[13][j] = coords->pv[11][j]; /* Nutations if on file */
			/*
			Calculate Earth and Moon from EMBary and geocentric Moon.
			*/
#if NewCode
			coords->pv[12][j] = coords->pv[2][j]; /* Move EMBary from Earth spot */
			coords->pv[2][j] -= coords->pv[9][j] / (1.0 + header->emrat); /* Earth */
			coords->pv[9][j] += coords->pv[2][j]; /* Moon (change geo->SS-centric) */
#else
			coords->pv[12][j] = coords->pv[2][j]; /* Move EMBary from Earth spot */
			coords->pv[2][j] -= coords->pv[9][j] / (1.0 + header->emrat); /* Earth */
			coords->pv[9][j] += coords->pv[2][j]; /* Moon (change geo->SS-centric) */
#endif
		}
		/*
		If we want heliocentric coordinates (not Solar System Barycenter
		coordinates), subtract Sun's position from all coordinates.
		*/
		if (!coords->bary) 
		{
			for (i=0; i<13; i++) 
			{
				if (i == coords->center) 
					i++; /* Skip over center object */
				if (i < 13) 
				{
					for (j=0; j<6; j++)
						coords->pv[i][j] -= coords->pv[coords->center][j];
				}
			}
			/*
			Set new center's positions and velocities to all zeroes
			*/
			coords->pv[coords->center][0] = coords->pv[coords->center][1] =
				coords->pv[coords->center][2] = coords->pv[coords->center][3] =
				coords->pv[coords->center][4] = coords->pv[coords->center][5] = 0.0;
		}
		
	}

	return(retval);
}


/*
ephcom_cheby() - interpolate at a point using Chebyshev coefficients
*/
int CEphemDE::ephcom_cheby(
						   double x,      /* Value of x over [-1,1] for Chebyshev interpolation */
						   double span,   /* Span in time of subinterval, for velocity */
						   double *y,     /* Chebyshev coefficients */
						   int ncoords,   /* Total number of coordinates to interpolate */
						   int ncoeffs,	/* Number of Chebyshev coefficients per coordinate */
						   double *pv,    /* Array to hold position in 1st half, velocity in 2nd */
						   bool bGetVelocity
						   ) 
{

	int i, j;
	double twox;
	int maxcoeffs = m_header.maxcheby;
	/* Position and velocity polynomial coefficients. */
	double * pc = &m_vecPos[0];
	double * vc = &m_vecVel[0];
	/*
	This need only be called once for each Julian Date,
	saving a lot of time initializing polynomial coefficients.
	*/
#if 0
	if (m_lastx != x) 
	{
		m_lastx = x;
		twox = x + x;   /* For Chebyshev recursion */
		/*
		Initialize position polynomial coefficients
		*/
		pc[0] = 1.0;    /* Chebyshev T[0](x) = 1 */
		pc[1] = x;      /* Chebyshev T[1](x) = x */
		// printf("pc[%d] = %20.13e\n", 0, pc[0]);
		// printf("pc[%d] = %20.13e\n", 1, pc[1]);
		{
			double *pc_0 = pc;
			double *pc_1 = pc + 1;
			double *pc_2 = pc + 2;
			for (i=2; i<maxcoeffs; i++, pc_0++, pc_1++, pc_2++) 
			{
				*pc_2 = twox * (*pc_1) - (*pc_0);
				/*
				Resolve bug with gcc generating -0.0 (also makes
				the smallest represented number equal to zero).
				*/
				if ((*pc_2) * (*pc_2) == 0.0) 
				{
					(*pc_2) = 0.0;
				}
			}
		}
	}
#else
	if (m_lastx != x)
	{
		m_lastx = x;
		twox = x + x;   
		/*
		Initialize position polynomial coefficients
		*/
		pc[0] = 1.0;    /* Chebyshev T[0](x) = 1 */
		pc[1] = x;      /* Chebyshev T[1](x) = x */
		// printf("pc[%d] = %20.13e\n", 0, pc[0]);
		// printf("pc[%d] = %20.13e\n", 1, pc[1]);
		for (i=2; i<maxcoeffs; i++) {
			pc[i] = twox * pc[i-1] - pc[i-2];

			if (pc[i]*pc[i] == 0.0) {
				pc[i] = 0.0;
			}
		}
		/*
		Initialize derivative polynomial coefficients
		*/
		vc[0] = 0.0;          /* d(1)/dx        = 0  */
		vc[1] = 1.0;          /* d(x)/dx        = 1  */
		vc[2] = twox + twox;  /* d(2x^2 - 1)/dx = 4x */
		for (i=3; i<maxcoeffs; i++) 
		{
			vc[i] = twox * vc[i-1] + pc[i-1] + pc[i-1] - vc[i-2];
			//    printf("vc[%2d] = 2(%g) * %g + 2 * %g - %g = %g\n",
			//           i, x, vc[i-1], pc[i-1], vc[i-2], vc[i]);
		}
	}
#endif

	if(bGetVelocity && m_lastxv != x)
	{
		m_lastxv = x;
		/*
		Initialize derivative polynomial coefficients
		*/
#if 0
		double *vc_0 = vc;
		double *vc_1 = vc + 1;
		double *vc_2 = vc + 2;

		double *pc_1 = pc + 1;

		*vc_0 = 0.0;          /* d(1)/dx        = 0  */
		*vc_1 = 1.0;          /* d(x)/dx        = 1  */
		*vc_2  = twox + twox; /* d(2x^2 - 1)/dx = 4x */
		for (i=3; i < maxcoeffs; i++, vc_0++, vc_1++, vc_2++, pc_1++) 
		{
			*vc_2 = twox * (*vc_1) + (*pc_1) + (*pc_1) - (*vc_0);
		}
		/*
		Initialize derivative polynomial coefficients
		*/
#else
// 		vc[0] = 0.0;          /* d(1)/dx        = 0  */
// 		vc[1] = 1.0;          /* d(x)/dx        = 1  */
// 		vc[2] = twox + twox;  /* d(2x^2 - 1)/dx = 4x */
// 		for (i=3; i < maxcoeffs; i++) 
// 		{
// 			vc[i] = twox * vc[i-1] + pc[i-1] + pc[i-1] - vc[i-2];
// 		}
#endif // 0
		//_ASSERT(0); // check optimization
	}
	/*
	Interpolate to get position for each component
	*/
#if 0
	for (i=0; i<ncoords; i++)  /* Once each for x, y, and z */
	{
		double *pPc = pc;
		double *pY = y + i*ncoeffs;
		double tmp = 0.0;
		for (j=ncoeffs-1; j >= 0; j--, ++pY, ++pPc) 
		{
			tmp += (*pPc) * (*pY);
		}
		pv[i] = tmp;
	}
#else
	for (i=0; i<ncoords; i++) { /* Once each for x, y, and z */
		pv[i] = 0.0;
		for (j=ncoeffs-1; j >= 0; j--) 
		{
			pv[i] += pc[j] * y[i*ncoeffs + j];
		}
	}
#endif
	/*
	Interpolate velocity (first derivative)
	*/
	if(bGetVelocity)
	{
#if 0
		for (i=0; i<ncoords; i++) 
		{
			double tmp = 0.0;
			double *pVc = vc;
			double *pY = y + i*ncoeffs;

			for (j=ncoeffs-1; j >= 0; j--, ++pY, ++pVc) 
			{
				tmp += (*pVc) * (*pY);
			}
			pv[ncoords + i] = tmp * 2.0 / span;
		}
#else
		for (i=0; i<ncoords; i++) {
			pv[ncoords + i] = 0.0;

			for (j=ncoeffs-1; j >= 0; j--) 
			{
				pv[ncoords + i] += vc[j] * y[i*ncoeffs + j];
			}
			pv[ncoords + i] *= 2.0 / span;
		}
#endif
	}
	else
	{
		for (i=0; i<ncoords; i++) 
		{
			pv[ncoords + i] = 0.0;
		}
	}
	return 0;
}

